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Abstract 

We use molecular dynamics simulations to study the 
formation of surface waves in vertically vibrated gran- 
ular material. We find that horizontal movements of 
particles, which are essential for the formation of the 
waves, consist of two distinct processes. First, the 
movements sharply increase while the particles are 
colliding with a bottom plate, where the duration of 
the collisions is very short compared to the period 
of the vibration. Next, the movements gradually de- 
crease between the collisions, during which the parti- 
cles move through the material. We also find that the 
horizontal velocity field after the collisions is strongly 
correlated to the surface profile before the collisions. 
PACS Number: 85.70.F: 47.20: 03.20. -|-i; 03.40.Kf 



Granular material under vibration has been a con- 
stant source of interesting phenomena among which 
surface waves are particularly interesting jl[ ||, ||, ^, 
^. Particles in a vertically vibrating box display 
a variety of steady state patterns — stripes, squares, 
hexagons, and localized excitations called "oscil- 



lons" , and even a richer set of transient patterns 
§, 0, |, g |l^, 0. The fact that such diverse 
patterns can exist in the seemingly simple system 
has attracted a lot of efforts on the study of the 
mechanism for the formation of the surface waves 
111, |ll,|ll0|7l|6l|8l|9|. 



It is generally accepted that horizontal movements 
of particles play an important role in the formation 
of the surface waves. The present theories can be 
roughly divided into two groups according to how 
such movements are generated and maintained. In 
the first group of theories, it is argued that horizontal 
movements are generated while particles are colliding 
with the bottom of a box |l|, ^ The duration 
of the collisions is very short compared to the pe- 
riod of the vibration. It is also argued that between 
the collisions with the bottom, horizontal movements 
gradually decrease while the particles move through 
the medium. Using these assumptions, many of the 
patterns observed in the experiments can be repro- 
duced. In the other group of theories, it is argued 
that horizontal movements depend only on contin- 
uum variables such as height and/or density fields 
p5| , p^ , p^ , p8| |. Since these fields change smoothly 
over time, changes in the movements are also grad- 
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ual. In particular, changes in the movements during 
the coUisions (with the bottom) are not particularly 
different from those during other phases of the vibra- 
tion. These theories can also reproduce many of the 
experimental patterns, typically by coupling density 
and height fields. 

The underlying assumptions of all of these theories 
sound plausible. Also, it is very difficult to deter- 
mine that which of the theories describes the experi- 
ments better. All of them reproduce many of the ex- 
perimental patterns, and have their own strong and 
weak points. Also, one should note that such theories 
cannot be validated just because they reproduce the 
observed patterns. Quantitative agreements (e.g., of 
dispersion relation) are necessary for the validation. 
In order to check the validity of the theories, it is 
thus necessary to obtain detailed information on the 
system, and compare it with the predictions of the 
theories. Unfortunately, such information is usually 
difficult to obtain by experiments. 

In this paper, we use molecular dynamics (MD) 
simulation method, which provides detailed informa- 
tion on the motion of individual particles as well as 
time averaged fields. We study the system in two 
dimensions. We focus on horizontal movements of 
particles, which are essential to the formation of the 
surface waves. We find that the time evolution of 
average horizontal speed U is made up of two sep- 
arate processes. First, there are sharp increases in 
U during the short periods that the pile is colliding 
with the bottom plate. The other process is grad- 
ual decay of U between such collisions, which results 
from interparticle collisions. The time evolution of 
U strongly supports the theories in the first group 
[ p2| , |l3|, . The present results do not imply that a 
continuum description is not possible for the system, 
it implies that the interpretation of the continuum 
variables has to be modified. 

We then study the processes in more detail. We 
find that the horizontal velocity field Vx{x) after col- 
lisions (with the bottom) shows a strong correlation 
with dxh{x) before the collisions, where h{x) is the 
center of mass field. Such correlation is assumed in 
the theories of 14 1. The second process can be 
characterized by temporal decay of U. We find that 
the decay time is very small when there is no sur- 



face wave, and it is comparable to the period of the 
vibration when surface waves are present. We also 
study the parameter dependence of the decay time 
and maximum horizontal speed. 

The simulations are done in two dimensions with 
disk shaped particles, using a form of interaction due 
to Cundall and Strack ||2^, . Particles interact only 
by contact, and the force between two such particles 
i and j is the following. Let the coordinate of the 
center of particle i (j) be Ri (Rj), and r = Ri — Rj- 
The normal component ^JLi of the force acting on 
particle i from particle j is 

F"_^i = kn{ai + a.j - \r\) - -fme{v-n), (1) 

where (a^) is the radius of particle i (j), n = f/r^ 
and V = dr/dt. Here, fc„ is the elastic constant, 7 
the friction coefficient, and me is the effective mass, 
rriimj / {mi+rrij). The shear component -FJ^^ is given 

by 



= -sign(5s) mm{k,\Ss\,^i\FJ'_ 



(2) 



where fi is the friction coefficient, ds the total shear 
displacement during a contact, and kg is the elastic 
constant of a virtual tangential spring. The shear 
force applies a torque to the particles, which then 
rotate. 

Particles can also interact with walls. The force 
and torque on particle i, in contact with a wall, are 
given by (|^) - (|^) with aj = and me = m^. Also, 
the system is in a vertical gravitational field g. The 
interaction parameters used in this study are fixed 
as follows, unless otherwise specified: kn = kg = 
10^,7 = 10'' and /i — 0.2. And, the timestep for 
integration is 5 x 10^^. In order to avoid artifacts 
of a monodisperse system (e.g., hexagonal packing), 
we choose the radius of the particles from a Gaussian 
distribution with the mean 0.1 and width 0.02. The 
density of the particles is 5. Throughout this paper, 
CGS units are implied. 

We put particles on a horizontal plate which os- 
cillates sinusoidally along the vertical direction with 
given amplitude A and frequency /. Let the width 
of the plate be W . We apply a periodic boundary 
condition in the horizontal direction. We start the 



2 



simulation by inserting the particles at random posi- 
tions above the plate. We let them fall by gravity and 
wait while they lose energy by collisions. We wait for 
10^ iterations for the particles to relax, and during 
this period we keep the plate fixed. The typical ve- 
locity at the end of the relaxation is of order 10~^. 
After the relaxation, we vibrate the plate, and start 
to take measurements. 

The coefficient of restitution between the particles 
Bpp, determined from the above interaction parame- 
ters, is 0.21, and the coefficient between the particles 
and the plate Cp^, is 8.0 x 10~^. The particles are thus 
strongly inelastic. We have studied the motion of a 
single particle for several values of A with / = 10, and 
find good agreements with the predictions by Mehta 
and Luck |22 . Also, the present model was shown 
to reproduce the dispersion relation from experiment 

il- 

We then study the center of mass motion of the 
particles. As the depth of the pile increases, its mo- 
tion can be different from that of a single particle 
p^ . For the parameters with which the motion of a 
single particle is periodic with period T = 1//, the 
motion of the pile can be subharmonic at large depth. 
Although studying the effect of the subharmonic mo- 
tion on the surface waves can be interesting, here we 
limit ourselves to the simpler case of no such motion. 
The minimum depth at which the subharmonic mo- 
tion occurs depends on the interaction parameters, 
and is an increasing function of k„. We thus use a 
rather large value of fc„(10^), so the subharnomic mo- 
tion is absent in all cases studied here. As a check, we 
study the motion of the particles in a narrow tube of 
W = 1 (five particle width) for several F with fixed /. 
Here, F is dimensionless peak acceleration AnAp /g. 
We find that the pile behaves like a single particle for 
at least 10 particle depth. In particular, a bifurcation 
of the motion occurs at F ~ 3.7, which then termi- 
nates at F ~ 4.4, which agree well with the results of 
a single inelastic particle [ p5| . 

We proceed to study surface waves. We choose 
VF = 16 and N = 800. Thus the system is, on av- 
erage, 80 particle wide and 10 particle deep. We fix 
the frequency / = 10, and study the waves for several 
values of F. Wc find that //2 waves start to appear 
for F ^ 2.5, then disappear when F becomes about 4. 



When F further increases, //4 waves start to appear 
for F ^ 5.5. The features of this "phase diagram" 
agree well with those of the experiments . 

Inspection of the motion of the particles shows that 
horizontal movements of the particles play a crucial 
role in maintaining the waves — an observation which 
was made in many previous studies on the problem 
(e.g, fill). We focus on the horizontal movements 
which we characterize by average horizontal speed 
U{t), defined as 



U{t) 



1 

N 



N 



(3) 



where u* is the horizontal velocity of particle i. Even 
if there are active horizontal movements, the average 
horizontal velocity can be small since the movements 
can occur in both positive and negative x directions. 
We thus use U{t) instead of the average horizontal ve- 
locity. We normalize U{t) with 2'KAf, the maximum 
velocity of the bottom plate. In Fig. 1, we show nor- 
malized U{t) for F = 2 and 3 for 10 vibration periods. 
Here, / = 10, VF = 16, and N = 800. The F = 3 
curve has been offset for clarity. 



Fig. 1 
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Figure 1: Dimensionless average horizontal speed 
U{t) for first 10 vibration cycles. Here, W = 16, 
N = 800, / = 10 and F = 2 and 3. The F = 3 curve 
has been offset for clarity. 

As shown in the figure, there clearly exist two sepa- 
rate processes: sharp increase in U within short time 
intervals and gradual decay of U during other phases 
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of the vibration. We also measure the time series of 
the total pressure on the plate p{t), which consists 
of sharp peaks occurring when the pile collides with 
the plate. We find that the locations of the peaks in 
U{t) coincide with those in p{t). Thus, the horizontal 
movements are being fed by collisions of the pile with 
the plate, and are being lost while the pile is not in 
contact with the plate. 

In order to maintain the surface waves, the par- 
ticles have to travel distance A within the period of 
the surface waves, where A is their wavelength. When 
U{t) decays slowly, we expect that the particles can 
travel long enough distance to maintain the waves. 
Indeed, we find that the surface waves are present 
only when U{t) decays slowly. To be more precise, 
the waves are observed when the decay time r, which 
will be defined later, is comparable to the period of 
the vibration. 

These observations strongly support the theories 
[ p2| , [T^ , which argue that the horizontal motion is 
being supplied by collisions with the bottom plate, 
and is being dissipated by interparticle collisions. 
The short time intervals during which the pile is col- 
liding with the plate play a crucial role in maintaining 
the surface waves. 



Fig. 2 




Figure 2: Time evolution of the dimensionless hori- 
zontal velocity field Vx{x, t). Here F = 3, / = 10, and 
other parameters are the same as in Fig. 1. The sud- 
den increases in Vx{x^t) occur when the pile collides 
with the plate. 

We then consider each process in more detail. In 



Fig. 2, we show time evolution of the horizontal ve- 
locity field Vx{x,t), which is defined as the average 
horizontal velocity of the particles whose centers are 
in [x, X -\- dx\. Here, F = 3, / = 10, and other param- 
eters are the same as in Fig. 1. We also normalize 
Vx{x,t) by 2Tr Af. In the figure, sharp changes in 
Vx{x,t) occur at certain values of t, which we find 
to be identical to the positions of the peaks in U{t). 
Thus, sudden increases in Vx{x,t) occur when the 
pile collides with the plate. The shape of Vx{x,t) af- 
ter such collisions does not change much until next 
collisions, but the overall magnitude of Vx{x,t) de- 
creases. 

The profile of Vx just after the collisions is impor- 
tant for the dynamics of the surface waves. Since the 
motion of the particles is deterministic, it is possible 
that the profile is determined from certain field just 
before the collisions. We check the dependency of Vx 
on several fields. Specifically, we calculate the corre- 
lation coefhcient r between Vx after the collisions and 
the following fields just before the collisions: x and y 
velocity fields Vx and Vy, their spatial derivative dxVx 
and dxVy, the number density m and its derivative 
dxin, and the center of mass field h and its derivative 
dxh. Here, m{x)dx {h{x)) is defined to be the total 
number (the center of mass) of the particles whose 
centers are in [a;,^ + dx\. Since the determination 
of the collision times becomes difficult for large F, we 
study only for F < 4. We average r over 10 collisions. 
We find that r is large (~ —0.8) for dxm, dxh and Vx, 
and there is essentially no correlation with the other 
fields. The large value of r for dxh suggests 

Vx{x,ta) oc -dxh{x,tb), (4) 

where ta and tb is the time just after and before the 
collisions, respectively. This very equation was as- 
sumed in the theories of Cerda et al. [|l2| and Eg- 
gers and Riecke jlj], where they proposed the form 
motivated from the flow of granular material on an 
inclined plane. The strong correlation occurs for all 
F studied here, including the F = 2 case where no 
surface wave is observed. 

We find that the spatial variation of the packing 
fraction is not significant, so m is roughly propor- 
tional to h. We thus expect that the correlation for 
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dxin is also large. What is strange is the strong cor- 
relation between Vx before and after the collisions, 
which states that the particles reverse the direction 
of their horizontal movements while the pile is collid- 
ing with the plate. We do not understand the origin 
of the reversal of the motion. 

We also measure the magnitude of the horizontal 
movements for several values of F. To be more spe- 
cific, we measure the peak values of U{t) like the 
ones shown in Fig. 1, and averaged them over all the 
collisions. One might expect that the average is pro- 
portional to the collision velocity of particles on the 
plate. Even though both display minimum around 
F — 4.6, we find that there is no strong correlation 
between them, and the average seems to show com- 
plicated dependence on F. Further study is necessary 
to quantify and understand the dependence. 

Next, we study the decay process, where horizon- 
tal movements of the particles gradually decrease be- 
tween the collisions with the plate. We quantify the 
process as follows. We start from a time series of U {t) 
like the ones in Fig. 1. We translate the positions of 
the peaks so that they all coincide. We then normal- 
ize the peak values of U (t) to unity, and average U {t) 
over all the peaks. The resulting quantity, defined as 
Ua{t), can be used to characterize the decay process. 
In Fig. 3, we show Ua{t) for F = 2,3 and 4, where 
other parameters remain unchanged. We then define 
decay time r as the "half-life" — the time at which 
Ua{t) becomes 1/2. The half-life for F = 2 is about 
0.0045, which is much smaller than the period of the 
vibration (T — 0.1). On the other hand, r is close to 
0.08 for F = 3 and 4. For all F wc have studied, we 
find that the surface waves are present if and only if 
r is comparable to the period of the vibration. The 
need for large r for the formation of the surface waves 
is evident as previously discussed. The fact that the 
dynamics consists of two separate processes does not 
seem to depend on small variations of interaction pa- 
rameters. 

The quantitative understanding of the decay pro- 
cess requires that one should be able to predict the 
dependence of r on parameters such as F. However, 
it seems that such dependences are rather complex. 
For example, we find that r does not depends on 
F in a simple fashion: r is not a monotonic func- 
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Figure 3: The time evolution of Ua{t) for F = 2, 3 and 
4. The "half-Hfe" r is small for F = 2, but becomes 
comparable to the period of the vibration for F = 3 
and 4. 

tion of F, and a small increase in F can result in 
a ten-fold increase or decrease in r. At the parti- 
cle level, the decay process occurs by the collisions 
between the particles, and the resulting changes in 
the horizontal movements. In order to make quanti- 
tative predictions, one thus has to understand both 
the collision frequency and the resulting momentum 
changes. Unfortunately, both of these quantities are 
poorly understood. More studies are again needed in 
order to understand the decay process. 

The time evolution of horizontal movements of the 
particles can also be studied using auto-correlation 
function c(t) defined as (vx{0)vx{t)), where the av- 
erage is taken over the particles. We find that c{t) 
also displays two basic processes: sharp change dur- 
ing the collisions with the plate and slow decay be- 
tween them. 

In sum, we show that horizontal movements of 
the particles, which are essential to the formation of 
the surface waves, consist of two separate processes: 
sharp increase of the movements by the collisions 
with the bottom plate, and slow decrease of the move- 
ments due to interparticle collisions, which strongly 
support the theories of ||l^, |l^, [l4| . We also find that 
the horizontal velocity field Vx {x) after the collisions 
with the plate is strongly correlated to dxh(x) just 
before the collisions. Here, we are mainly interested 
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in qualitative features of the mechanism for the for- 
mation of the surface waves. Their quantitative un- 
derstanding, such as the parameter dependence of r, 
will be subject of future works. 

After the present work has been completed, we be- 
came aware of the work of Kim et al p7[ | , where sim- 
ilar results were obtained. 

I thank H. K. Pak, S. Kim, K. Kim and S.-O. Jeong 
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by the Department of Energy under grant DE-FG02- 
93-ER14327, Korea Science and Engineering Founda- 
tion through the Brain-Pool program, and SNU-CTP. 
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